Evolution of inspiral orbits around a Schwarzschild black hole 
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We present results from calculations of the orbital evolution in eccentric binaries of nonrotating black holes 
with extreme mass-ratios. Our inspiral model is based on the method of osculating geodesies, and is the first to 
incorporate the full gravitational self-force (GSF) effect, including conservative corrections. The GSF informa- 
tion is encapsulated in an analytic interpolation formula based on numerical GSF data for over a thousand sample 
geodesic orbits. We assess the importance of including conservative GSF corrections in waveform models for 
gravitational-wave searches. 
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Introduction. The relativistic two-body problem for bina- 
ries of greatly different masses is of both theoretical impor- 
tance and astrophysical relevance. When the mass ratio is 
extreme, deviation from geodesic motion can be described 
in terms of an effective gravitational self-force (GSF) arising 
from the interaction of the small object with its own space- 
time perturbation. The fundamental problem of regularizing 
the gravitational self-interaction in curved spacetime has been 
studied extensively since the late 1990s, and is now rigorously 
solved at the first post-geodesic order (i.e., at linear order in 
the small mass ratio rj) 0J31- This theoretical advance has 
been strongly motivated by the exciting prospects of observ- 
ing gravitational waves from inspiralling black hole binaries 
of small mass ratios. Systems with rj = 10 2 — 10 3 (IMRIs: 
intermediate-mass-ratio inspirals) could be detected by Ad- 
vanced LIGO [5 1, and a low -frequency detector in space based 
on the LISA design will observe many inspiralling systems 
with rj — 10~ 4 -10~ 7 (EMRIs: extreme-mass-ratio inspirals) 
[6 1. The latter, involving a compact object captured by a mas- 
sive black hole, are of key importance in gravitational-wave 
astronomy due to their unique utility as probes of strong-field 
gravity Q IS) . To interpret the information encoded in the 
I/EMRI signals it is crucial to have at hand an accurate model 
of the orbital evolution driven by the GSF. 

The last decade has seen a concentrated effort to develop 
computational tools for the GSF in black hole spacetimes 
||9| . This culminated in 2010 [ 1 1 with the introduction of 
a code that returns the GSF along any specified (fixed) bound 
geodesic orbit of the Schwarzschild geometry. Results from 
this code have already been used to quantify aspects of the 
conservative post-geodesic dynamics ifTTI . and today they 
provide a unique strong-field benchmark for post-Newtonian 
calculations fl2l . nonlinear numerical simulations ifTSl [T4l . 
and Effective One Body theory lfl5HT71 . However, with the 
I/EMRI problem in mind, it remains an important task to 
translate the GSF information into inspiral trajectories and 
gravitational waveforms. This has not been attempted so far. 

Here we report an important milestone in the GSF pro- 
gramme: an algorithm and a working code for computing 
inspiral orbits on a Schwarzschild background, incorporating 
the full GSF information. The full GSF has a dissipative piece, 



responsible for the orbital decay, but also a conservative com- 
ponent, which, e.g., modifies the rate of periastron precession. 
In the radiative approximation (RA) one ignores the conserva- 
tive GSF effect and considers only the secular, time-averaged 
part of the dissipative dynamics. The latter can be computed 
using global energy-momentum balance considerations, with- 
out resorting to the local GSF [18|. The RA can bring con- 
siderable computational saving, so it is important to assess its 
efficacy, which we do here reliably for the first time. 

There are two approaches to self-forced evolution. In the 
systematic approach one solves the perturbation equations and 
the self-forced equations of motion as a coupled set, in a self- 
consistent manner. This entails incorporating back-reaction 
corrections in the GSF code itself — a technically challenging 
task yet to be attempted. The second approach invokes the 
traditional method of osculating orbits (also known in New- 
tonian celestial mechanics as the method of variation of con- 
stants). In this approach the inspiral orbit is reconstructed as a 
smooth sequence of geodesies, each lying tangent to the orbit 
at a particular moment. This amounts to modelling the true or- 
bit as an evolving geodesic with dynamical orbital elements. 
Equations governing the forced evolution of the latter in the 
Schwarzschild case (with an arbitrary forcing agent) were ob- 
tained by Pound and Poisson lfT9l . and Gair et al. generalized 
the formalism to Kerr 11201 . We adopt here the formalism of 
|l9l , and implement it with actual GSF data for the first time. 

Throughout this letter we set G = c = 1 and use metric sig- 
nature ( — h++) and Schwarzschild coordinates {t,r,6,q>}. 
M denotes the mass of the background Schwarzschild geome- 
try, and jj. is the mass of the inspiralling object (so /J./M = rj). 

Osculating geodesies. Bound geodesies of the 
Schwarzschild geometry can be parametrized by their 
semilatus rectum pM and eccentricity <?, defined via 
r± = pM I (1 =F e), where r = r + and r = r_ are the apastron 
and periastron radii, respectively. The geodesic motion of a 
test particle is described by 
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where %(t) is a monotonically increasing parameter along the 
orbit, obtained by inverting t(%) — §~Sdt j d%')d%' with 
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Here v = %— Xo, and without loss of generality we assumed 
the motion takes place in the equatorial plane (9 — 7C/2), and 
t (Xo) — f(Zo) = ( f = is a periastron passage with <p = 0). 

In the osculating geodesies approach, the inspiral motion 
of a mass particle under the effect of the GSF is described 
by r = r g {f,p{t),e{t),Xo{t)) and <p = <p g (f,p(t),e(t),Xo{t)), 
where pit), e(t), Xo{t) are osculating elements. The principal 
elements p and e determine the "shape" of the orbit; the posi- 
tional element Xo describes the orientation of the major axis. 
Both principal and positional elements evolve secularly under 
the GSF effect, but while the secular evolution of p and e is 
dissipative, that of Xo is conservative (it describes the preces- 
sion effect of the GSF). Both principal and positional elements 
also exhibit quasi-periodic oscillations. 

Given the GSF components F a (°c T7 2 ), the osculating ele- 
ments evolve according to Hl9l 
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where an overdot denotes d/dt, F a = jti l F a is the self- 
acceleration, /o = (p — 2 — 2ecosv)(p — 3 — e 2 )[(p — 2) 2 — 
4e 2 ]- 1 /2[Q,_ 6 )2_ 4e 2]-i i ft = ( / ,_6-2ecosv) 1 / 2 , f 2 = 

(l+ecosv)" 2 ,/ 3 =/ 2 ecosv + 2(/?-3) and J3 = p - 6 - 2e 2 . 
In our implementation, the GSF data will be given in the 
form F a = F a (x — Xo', P, e ), evaluated along geodesies with 
fixed p,e,Xo- With these data at hand, Eqs. |4]l-((6]l [with Eq. 
Q] form a closed set of ordinary differential equations for 
{p(t),e(t),Xo{t )}■ We will solve this set with the initial condi- 
tions {p(0),e(0),^o(0)} = {po,e ,Q} for some po,e () . The in- 
spiral trajectory will then be described by Eqs. ([TJ1 and (|2| with 
p,e,Xo replaced by the corresponding osculating elements. 

GSF interpolation model. Existing codes do not return 
the true GSF along the evolving orbit, but an approximation 
thereof computed along fixed geodesies. The resulting error is 
very small in the adiabatic regime where the evolution occurs 
on a timescale much longer than the orbital period T. It can be 
shown |21| that the adiabaticity condition a = (\p/p\T) <C 1 
(where (•) denotes an average over time T) is met so long as 
£ S> 17 1//2 , where e = p — 6 — 2e is a measure of the proximity 
to the innermost stable orbit (ISO). Thus, for EMRI-relevant 
tj values, the evolution is adiabatic until very close to the ISO. 
Beyond that point, our GSF model may cease to be useful. 



GSF codes return F a (x) for given p,e,Xo- To express 
this information in a workable form we devise accurate an- 
alytic fits to numerical data obtained using two independent 
codes: the original, time-domain code of Ref. [ 10], and a new 
code based on a frequency-domain treatment of the Lorenz- 
gauge perturbation equations l22l l23l . Both codes take as 
input the geodesic parameters p,e, and return (separately) 
the dissipative and conservative pieces of the GSF along the 
geodesic, F^ ss (X',P,e) and F c ® ns (x;p,e) respectively. The 
new, frequency-domain, algorithm offers significant compu- 
tational saving, particularly at low eccentricity. This is a cru- 
cial improvement, since GSF calculations are extremely com- 
putationally intensive. We used our frequency-domain code 
to compute the GSF for a dense sample of p,e values in the 
range < e < 0.2 and 6 + 2e < p < 12. We tested a subset of 
the results using our time-domain code. The results described 
below are based on a sample of 1 100 geodesies, for which the 
GSF has been computed with fractional accuracy < 10 -4 . 

To devise an interpolation formula for the numerical data, 
we observe that the GSF is a periodic function of v = % — Xo 
along the geodesic, with F^,F c r ons even in v, and F d r iss ,F c ^ ns 
odd in v (9). This suggests the Fourier-like representation 
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(i = 1,...,4), where J?- = J?- 2 ^,^^^^^} 
and osCj(nv) = {cosnv, sinnv, sinnv, cosnv}. We have em- 
ployed here a simple power series model for the p,e depen- 
dence of the Fourier coefficients. For the leading l/p power 
we take kj = {2, 2,4, -y}, consistent with the known behav- 
ior of the various ^,'s at large p. Each Fourier «-mode of ^ 
admits a power series in e 2 starting at e n . The dimensionless 
numerical coefficients a, n ^ in Eq. (|7jl are to be determined 
by fitting to numerical GSF data, with the summation cutoffs 
Hi,ji,ki to be chosen empirically. 

We used a standard least-squares algorithm to fit the inter- 
polation formula (|7]i to the numerical data over the range of 
p,e-values indicated above. For our illustrative computation 
we sought a fractional accuracy < 10~ 3 in each component J?; 
[i.e., we demanded that Eq. (|7]) reproduced all available data 
to within that accuracy]. We found empirically that this can 
be achieved with n,- = 6, ji = 2 and k{ = 9 for each i. Thus the 
procedure fits 7 x 3 x 10 = 210 parameters a,„^. using 1100 
data points for each i. For lack of space we do not give here 
the values of the best-fit model parameters flj„ ;J t, but we have 
made them available online on a dedicated website [24] as part 
of an open-source "fast GSF calculator". The package con- 
tains a script for computing the GSF quickly based on Eq. (|7]) 
and a database of a;„» coefficients. We intend to update the 
database regularly as more GSF data (of improved accuracy 
and greater extent in the p,e space) become available. 

Sample results. Figures[T|and|2]display results from a full- 
GSF inspiral with rj = 10~ 5 , starting at (po,e ) = (12,0.2) 
(and taking M = 1O 6 M for concreteness). The orbit decays 
adiabatically (see the lower inset in Fig. [2} and circularizes 
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gradually, until very close to the ISO where the eccentricity 
begins to increase — a phenomenon already described in Ref. 
11251 . The entire inspiral, from po = 12 to the onset of plunge, 
lasts ~ 1443 x M/(1O 6 M ) days, during which the orbit com- 
pletes 75,550 periastron passages. Note the periastron phase 
Xq shifts secularly in a retrograde sense (in our example, by 
~ 9 radians over the entire inspiral). This represents a GSF- 
induced decrease in the rate of periastron advance (cf. [11 J). 
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FIG. 1. Sample full-GSF inspiral orbit with fi = 10M @ and M = I0''W. . 
starting at (pg,eo) = (12,0.2). We plot the orbit in the plane of x = 
(r/M) cos (p and y = (r/M) sin <p, showing 4 episodes during the inspiral: the 
onset of inspiral (~ 1443 days to plunge; top left), 500 days to plunge (top 
right), 75 days to plunge (bottom left), and the last hour of inspiral (bot- 
tom right). The motion is counter-clockwise, and each of the tracks shows 
one hour of inspiral. The central black hole (but not the orbiter) is drawn to 
scale. Periastron passages are indicated along with their sequential number, 
counting from the initial periastron (W = 0'). In the last snapshot the orbit 
completes ~ 6.7 revolutions in <p between the two periastra shown. 

Radiative approximation. To explore the long-term effect 
of the GSF's conservative piece, let us construct an RA model 
by setting F c " ns = in the evolution equations (|4|-((6|, and 
additionally replacing the expressions on the right-hand side 
with their corresponding f-averages over an entire radial pe- 
riod of the instantaneous osculating geodesic. We ask how 
well this RA model can capture the full-GSF dynamics. 

As a reference for comparison let us consider the accumu- 
lated azimuthal phase <p(f). We denote by (Pra/mi the values 
corresponding to the RA/full GSF models, and aim to inspect 
how the phase difference A^Jra = <Pra — <Pfuii builds up over 
time. To define A^Jra unambiguously we must map cautiously 
between the initial parameters of the RA and full-GSF models, 
noting the 0{r\ ) gauge ambiguity in the values of the param- 
eters p,e. A mapping based on "same po,eQ values" would 
result in the RA and full orbits possessing different initial fre- 
quencies, because the conservative piece of the GSF, which is 
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FIG. 2. Evolution of the osculating elements in the sample case of Fig. [T] 
We show the eccentricity e [lighter (red) line, left axis] and periastron phase 
Xq [darker (blue) line, right axis] as functions of semi-latus rectum p, as the 
binary inspirals from po = 12 down to the ISO (dashed curve). Marks along 
the curves count down (from right to left) 500 days, 100 days, 10 days, 1 day 
and 1 hour to the onset of plunge. Note the orbit initially circularizes, but 
upon approaching the ISO the eccentricity begins to increase. Note also the 
phase Xo decreases monotonically, implying that the conservative GSF acts to 
reduce the rate of relativistic precession. The upper inset is an enlargement of 
the near-ISO region; the manifest oscillatory behavior is due to the variation 
of the GSF with the radial phase. The lower inset shows the magnitude of 
the adiabaticity parameter a = (\p/p\T) vs. the ISO distance £ = p — 6 — 2e, 
confirming that the evolution is strongly adiabatic until very near the ISO. 



accounted for in the full model but not in RA, shifts the fre- 
quencies by an amount of 0(r\). This, in turn, would result in 
a rapid linear-in-time growth of AcjJra. 

To eliminate this spurious effect we instead match the fre- 
quencies of the initial osculating geodesies, using knowledge 
of F c " ns . To achieve this in practice we apply the following 
procedure. (1) Choose po,eo for the full-GSF orbit (taking 
%q = <p = at t = 0). (2) Compute the azimuthal and ra- 
dial frequencies of the orbit at t = through 0(rj), includ- 
ing F c " ns -induced corrections. (3) Find the p,e values of a 
geodesic whose frequencies are those found in step 2. (4) Use 
these p,e as initial values for the RA evolution (starting again 
with Xo = <P = at t = 0). We explain this procedure in more 
detail in a follow-up work ||26ll , which further explores the 
performance of the RA model. Our procedure matches the 
initial frequencies of the full and RA orbits. This is physically 
motivated because the frequencies (unlike p,e) are invariant 
characteristics of the orbit. 

Figure[3]shows A<Pra(?) for our sample orbit with tj = 1CT 5 
and (po,eo) = (12,0.2). On the lower horizontal axis we ex- 
press t in units of the radiation-reaction timescale ?rr = T c / T], 
where as a characteristic orbital period we take the ^-period 
of the innermost stable circular orbit, T c = 2%(y'l 1 M. As ex- 
pected, A<Pra grows secularly in proportion to (?/?rr) 2 (with 
oscillations reflecting the mismatch in radial phase between 
the RA and full-GSF orbits). This secular growth is attributed 
to conservative corrections to the rate-of-change of the az- 
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imuthal frequency, which are 0{r\ 2 ). The phase difference 
AcjjRA remains small for quite long, becoming significant only 
on a timescale t ~ ?rr. For reference, we also show in Fig. [3] 
the phase difference without adjusting the initial frequencies, 
i.e., using the same po, eo values for both models. In this case 
A<Pra x t, and (jJra quickly drifts away with respect to <p m n. 

t (days) 
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FIG. 3. Effect of conservative GSF corrections on the long-term phase evo- 
lution. Plotted is the accumulated phase difference A(J)r a = <j>ra — <Pf u ii (RA: 
Radiative Approximation) for the sample orbit shown in Fig.^ In the lower 
(blue) curve we have matched the initial frequencies of the RA and full orbits 
[correcting for the initial 0(r\) frequency shift due to the conservative piece 
of the GSF]. The upper (red) curve shows, for reference, the phase difference 
APra when no such adjustment is made. 

The RA model appears to capture the full-GSF phase evo- 
lution rather well over an extended portion of the inspi- 
ral (allowing for a suitable adjustment of the initial param- 
eters). Our results confirm the expectation that RA-based 
waveform templates could be implemented usefully in semi- 
coherent matched filtering searches for gravitational waves 
from I/EMRIs l27ij29l To obtain a fully phase-coherent the- 
oretical model of the evolution beyond the radiation-reaction 
timescale requires the conservative GSF, but it also requires 
the (as yet unknown) second-order piece of the dissipative 
GSF. Secular geodesic effects associated with the spin of the 
small object may also be important at this order. 

Concluding remarks. We reported here the development of 
a computational framework for calculating fully self-forced 
inspiral orbits for I/EMRI applications. This framework is un- 
der continuing development, and several important improve- 
ments must be made before we can compute fully coherent 
waveforms for astrophysical I/EMRIs. First, more accurate 
GSF data must be obtained and implemented to inform a 
more accurate interpolation formula; a fractional accuracy of 
< 0(r\) in the GSF is desired to ensure that the GSF model 
error has a negligible long-term effect. Such an accuracy stan- 
dard is achievable in principle using current codes El l23l . 
but advanced computational techniques now being developed 
will allow crucial saving in computational cost. Second, our 
GSF model must include 2nd-order dissipative corrections. 
Work to formulate and compute such corrections is under way. 



Third, the model must be extended to Kerr geometry. Tech- 
niques to compute the GSF in Kerr are under active devel- 
opment QUI [3D . Once Kerr GSF data are at hand, an inter- 
polation model akin to |7]i will need to be devised and im- 
plemented. Finally, it remains to quantify the error coming 
from calculating the GSF along fixed geodesies (and not the 
true evolving orbit). This must await the completion of a fully 
self-consistent evolution code, also under development. 

It would be instructive to compare our inspiral orbits with 
results from fully nonlinear numerical simulations once these 
become available for small r\. The state of the art is a short 
simulation for rj = 1:100 [32], covering the last few orbits of 
inspiral. A much longer simulation would be required to allow 
a meaningful comparison with our GSF results. 
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